# This file produces figures in the text and supplementary information for
#   Goehring and Lowande's Public Responses to Unilateral Policymaking. 

# See the README for more information. 

# load required packages and data
require(arm)
require(lubridate)
require(gridExtra)
require(ordinal)
require(stargazer)
require(kableExtra)
require(mfx)
require(xtable)
require(slider)
require(tidyverse)
require(broom)
require(irr)
require(MASS)

load('input-data/complete-data-v3.Rda') 

########################
## figure 2, figure a3 #
########################

# TOPIC APPROVAL
t1=clm(topic.approve.w1~topic.w1+president.w1+condition.w1,data=dt,family='logit')
t2=clm(topic.approve.w1~topic.w1+president.w1+partisanship+sex+age+income+ethnicity+condition.w1,data=dt,family='logit')
t3=clm(topic.approve.w1~topic.w1+president.w1+partisanship+sex+age+income+education+ethnicity+condition.w1,data=dt,family='logit')

# PRESIDENTIAL APPROVAL
p1=clm(pres.approve.w1~topic.w1+president.w1+condition.w1,data=dt,family='logit')
p2=clm(pres.approve.w1~topic.w1+president.w1+partisanship+sex+age+income+ethnicity+condition.w1,data=dt,family='logit')
p3=clm(pres.approve.w1~topic.w1+president.w1+partisanship+sex+age+income+education+ethnicity+condition.w1,data=dt,family='logit')

# TRUMP VOTE
v1=glm(trumpvote.w1~topic.w1+president.w1+condition.w1,data=dt,family='binomial')
v2=glm(trumpvote.w1~topic.w1+president.w1+partisanship+sex+age+income+ethnicity+condition.w1,data=dt,family='binomial')
v3=glm(trumpvote.w1~topic.w1+president.w1+partisanship+sex+age+income+education+ethnicity+condition.w1,data=dt,family='binomial')

# WORKS WITH CONGRESS
w1=clm(manip.cong.w1~topic.w1+president.w1+condition.w1,data=dt,family='logit')
w2=clm(manip.cong.w1~topic.w1+president.w1+partisanship+sex+age+income+ethnicity+condition.w1,data=dt,family='logit')
w3=clm(manip.cong.w1~topic.w1+president.w1+partisanship+sex+age+income+education+ethnicity+condition.w1,data=dt,family='logit')

# GETS THINGS DONE
g1=clm(manip.done.w1~topic.w1+president.w1+condition.w1,data=dt,family='logit')
g2=clm(manip.done.w1~topic.w1+president.w1+partisanship+sex+age+income+ethnicity+condition.w1,data=dt,family='logit')
g3=clm(manip.done.w1~topic.w1+president.w1+partisanship+sex+age+income+education+ethnicity+condition.w1,data=dt,family='logit')

# CARES ABOUT THE CONSTITUTION AND RULE OF LAW
l1=clm(manip.law.w1~topic.w1+president.w1+condition.w1,data=dt,family='logit')
l2=clm(manip.law.w1~topic.w1+president.w1+partisanship+sex+age+income+ethnicity+condition.w1,data=dt,family='logit')
l3=clm(manip.law.w1~topic.w1+president.w1+partisanship+sex+age+income+education+ethnicity+condition.w1,data=dt,family='logit')

# Estimates for main figure

# re-run as logits
t3.bin=glm(topic.approve.w1.bin~topic.w1+president.w1+partisanship+sex+age+income+education+ethnicity+condition.w1,data=dt,family='binomial')
p3.bin=glm(pres.approve.w1.bin~topic.w1+president.w1+partisanship+sex+age+income+education+ethnicity+condition.w1,data=dt,family='binomial')
c3.bin=glm(manip.cong.w1.bin~topic.w1+president.w1+partisanship+sex+age+income+education+ethnicity+condition.w1,data=dt,family='binomial')
g3.bin=glm(manip.done.w1.bin~topic.w1+president.w1+partisanship+sex+age+income+education+ethnicity+condition.w1,data=dt,family='binomial')
l3.bin=glm(manip.law.w1.bin~topic.w1+president.w1+partisanship+sex+age+income+education+ethnicity+condition.w1,data=dt,family='binomial')

outcomes=matrix(data=NA, ncol=7,nrow=18)

public_comparisons=function(model,data,simulations=1000,seed=6658) {
  set.seed(seed)
  betas=slot(sim(model,simulations),'coef')
  covariates=names(model$model)
  results=matrix(data=NA, ncol=7,nrow=3)
  
  # add columns for factor coefs
  D=data[complete.cases(data[,covariates]),]
  D=cbind(D,model.matrix(~topic.w1-1,data=D))
  D=cbind(D,model.matrix(~president.w1-1,data=D))
  if ('partisanship' %in% covariates) {D=cbind(D,model.matrix(~partisanship-1,data=D))}
  if ('sex' %in% covariates) {D=cbind(D,model.matrix(~sex-1,data=D))}
  if ('income' %in% covariates) {D=cbind(D,model.matrix(~income-1,data=D))}
  if ('education' %in% covariates) {D=cbind(D,model.matrix(~education-1,data=D))}
  if ('ethnicity' %in% covariates) {D=cbind(D,model.matrix(~ethnicity-1,data=D))}
  D=cbind(D,model.matrix(~condition.w1-1,data=D))
  
  # intercept + covariates 
  controls=cbind(1,D$`topic.w1contractor pay`,D$`topic.w1farm payments`,D$`topic.w1foreign aid and abortion`,D$`topic.w1foreign workers`,D$`topic.w1gun violence research`,D$topic.w1LGBT,D$`topic.w1policing and military surplus`,D$`topic.w1public lands`,D$`topic.w1Russian sanctions`,D$`topic.w1student loans`,D$topic.w1trade,D$`topic.w1water rules`,D$`topic.w1weapon sales to Saudi Arabia`,D$topic.w1wildlife,D$president.w1Trump)
  if ('partisanship' %in% covariates) {controls=cbind(controls,D$partisanshipDemocrat,D$partisanshipRepublican)}
  if ('sex' %in% covariates) {controls=cbind(controls,D$sexFemale)}
  if ('age' %in% covariates) {controls=cbind(controls,D$age)}
  if ('income' %in% covariates) {controls=cbind(controls,D$`income2: $25,000 to $50,000`,D$`income3: $50,001 to $75,000`,D$`income4: $75,001 to $100,000`,D$`income5: More than $100,001`)}
  if ('education' %in% covariates) {controls=cbind(controls,D$`educationSome College or Vocational`,D$`educationB.A. or B.S.`,D$`educationPost-graduate or Higher`)}
  if ('ethnicity' %in% covariates) {controls=cbind(controls,D$ethnicityBlack,D$`ethnicityNative American`,D$`ethnicityOther/Decline to State`,D$ethnicityWhite)}
  
  # new data to simulate over
  d.position=cbind(controls,0,0)
  d.order=cbind(controls,0,1)
  d.congress=cbind(controls,1,0)
  
  # simulations
  p.pos=apply((1/(1+exp(-as.matrix(d.position)%*% t(betas)))),1,as.vector) 
  p_pos=apply(p.pos,1,mean)
  p.ord=apply((1/(1+exp(-as.matrix(d.order)%*% t(betas)))),1,as.vector) 
  p_ord=apply(p.ord,1,mean)
  p.cong=apply((1/(1+exp(-as.matrix(d.congress)%*% t(betas)))),1,as.vector) 
  p_cong=apply(p.cong,1,mean)
  
  # store
  results[1,1]<- quantile(p_ord-p_pos,.0083)
  results[1,2]<- quantile(p_ord-p_pos,.025)
  results[1,3]<- mean(p_ord-p_pos)
  results[1,4]<- quantile(p_ord-p_pos,.975)
  results[1,5]<- quantile(p_ord-p_pos,.9916)
  results[1,6]<- 'Executive Order'
  results[1,7]<- covariates[1]
  results[2,1]<- quantile(p_cong-p_pos,.0083)
  results[2,2]<- quantile(p_cong-p_pos,.025)
  results[2,3]<- mean(p_cong-p_pos)
  results[2,4]<- quantile(p_cong-p_pos,.975)
  results[2,5]<- quantile(p_cong-p_pos,.9916) 
  results[2,6]<- 'Congress'
  results[2,7]<- covariates[1]
  results[3,1]<- quantile(p_ord-p_cong,.0083)
  results[3,2]<- quantile(p_ord-p_cong,.025)
  results[3,3]<- mean(p_ord-p_cong)
  results[3,4]<- quantile(p_ord-p_cong,.975)
  results[3,5]<- quantile(p_ord-p_cong,.9916)
  results[3,6]<- 'Executive Order - Congress'
  results[3,7]<- covariates[1]
  
  return(results)
}

outcomes=data.frame(lb.MC=numeric(18),lb.95=numeric(18),me=numeric(18),ub.95=numeric(18),ub.MC=numeric(18),comparison=character(18),dv=character(18))
outcomes[1:3,1:7]=public_comparisons(t3.bin,dt)
outcomes[4:6,1:7]=public_comparisons(p3.bin,dt)
outcomes[7:9,1:7]=public_comparisons(v3,dt)
outcomes[10:12,1:7]=public_comparisons(c3.bin,dt)
outcomes[13:15,1:7]=public_comparisons(g3.bin,dt)
outcomes[16:18,1:7]=public_comparisons(l3.bin,dt)
outcomes_overall=outcomes

# hypothesis 1A: break down co-partisanship

# TOPIC APPROVAL
t1a=clm(topic.approve.w1~topic.w1+president.w1+copartisanship*condition.w1,data=dt,family='logit')
t2a=clm(topic.approve.w1~topic.w1+president.w1+sex+age+income+ethnicity+copartisanship*condition.w1,data=dt,family='logit')
t3a=clm(topic.approve.w1~topic.w1+president.w1+sex+age+income+education+ethnicity+copartisanship*condition.w1,data=dt,family='logit')

# PRESIDENTIAL APPROVAL
p1a=clm(pres.approve.w1~topic.w1+president.w1+copartisanship*condition.w1,data=dt,family='logit')
p2a=clm(pres.approve.w1~topic.w1+president.w1+sex+age+income+ethnicity+copartisanship*condition.w1,data=dt,family='logit')
p3a=clm(pres.approve.w1~topic.w1+president.w1+sex+age+income+education+ethnicity+copartisanship*condition.w1,data=dt,family='logit')

# TRUMP VOTE
v1a=glm(trumpvote.w1~topic.w1+president.w1+copartisanship*condition.w1,data=dt,family='binomial')
v2a=glm(trumpvote.w1~topic.w1+president.w1+sex+age+income+ethnicity+copartisanship*condition.w1,data=dt,family='binomial')
v3a=glm(trumpvote.w1~topic.w1+president.w1+sex+age+income+education+ethnicity+copartisanship*condition.w1,data=dt,family='binomial')

# WORKS WITH CONGRESS
w1a=clm(manip.cong.w1~topic.w1+president.w1+copartisanship*condition.w1,data=dt,family='logit')
w2a=clm(manip.cong.w1~topic.w1+president.w1+sex+age+income+ethnicity+copartisanship*condition.w1,data=dt,family='logit')
w3a=clm(manip.cong.w1~topic.w1+president.w1+sex+age+income+education+ethnicity+copartisanship*condition.w1,data=dt,family='logit')

# GETS THINGS DONE
g1a=clm(manip.done.w1~topic.w1+president.w1+copartisanship*condition.w1,data=dt,family='logit')
g2a=clm(manip.done.w1~topic.w1+president.w1+sex+age+income+ethnicity+copartisanship*condition.w1,data=dt,family='logit')
g3a=clm(manip.done.w1~topic.w1+president.w1+sex+age+income+education+ethnicity+copartisanship*condition.w1,data=dt,family='logit')

# CARES ABOUT THE CONSTITUTION AND RULE OF LAW
l1a=clm(manip.law.w1~topic.w1+president.w1+copartisanship*condition.w1,data=dt,family='logit')
l2a=clm(manip.law.w1~topic.w1+president.w1+sex+age+income+ethnicity+copartisanship*condition.w1,data=dt,family='logit')
l3a=clm(manip.law.w1~topic.w1+president.w1+sex+age+income+education+ethnicity+copartisanship*condition.w1,data=dt,family='logit')

# re-run as logits
t3a.bin=glm(topic.approve.w1.bin~topic.w1+president.w1+sex+age+income+education+ethnicity+copartisanship*condition.w1,data=dt,family='binomial')
p3a.bin=glm(pres.approve.w1.bin~topic.w1+president.w1+sex+age+income+education+ethnicity+copartisanship*condition.w1,data=dt,family='binomial')
v3a=glm(trumpvote.w1~topic.w1+president.w1+sex+age+income+education+ethnicity+copartisanship*condition.w1,data=dt,family='binomial')
c3a.bin=glm(manip.cong.w1.bin~topic.w1+president.w1+sex+age+income+education+ethnicity+copartisanship*condition.w1,data=dt,family='binomial')
g3a.bin=glm(manip.done.w1.bin~topic.w1+president.w1+sex+age+income+education+ethnicity+copartisanship*condition.w1,data=dt,family='binomial')
l3a.bin=glm(manip.law.w1.bin~topic.w1+president.w1+sex+age+income+education+ethnicity+copartisanship*condition.w1,data=dt,family='binomial')

public_comparisons_party=function(model,data,simulations=1000,seed=6658) {
  set.seed(seed)
  betas=slot(sim(model,simulations),'coef')
  covariates=names(model$model)
  results=matrix(data=NA, ncol=8,nrow=9)
  
  # add columns for factor coefs
  D=data[complete.cases(data[,covariates]),]
  D=cbind(D,model.matrix(~topic.w1-1,data=D))
  D=cbind(D,model.matrix(~president.w1-1,data=D))
  if ('copartisanship' %in% covariates) {D=cbind(D,model.matrix(~copartisanship-1,data=D))}
  if ('sex' %in% covariates) {D=cbind(D,model.matrix(~sex-1,data=D))}
  if ('income' %in% covariates) {D=cbind(D,model.matrix(~income-1,data=D))}
  if ('education' %in% covariates) {D=cbind(D,model.matrix(~education-1,data=D))}
  if ('ethnicity' %in% covariates) {D=cbind(D,model.matrix(~ethnicity-1,data=D))}
  D=cbind(D,model.matrix(~condition.w1-1,data=D))
  
  # intercept + covariates 
  controls=cbind(1,D$`topic.w1contractor pay`,D$`topic.w1farm payments`,D$`topic.w1foreign aid and abortion`,D$`topic.w1foreign workers`,D$`topic.w1gun violence research`,D$topic.w1LGBT,D$`topic.w1policing and military surplus`,D$`topic.w1public lands`,D$`topic.w1Russian sanctions`,D$`topic.w1student loans`,D$topic.w1trade,D$`topic.w1water rules`,D$`topic.w1weapon sales to Saudi Arabia`,D$topic.w1wildlife,D$president.w1Trump)
  if ('sex' %in% covariates) {controls=cbind(controls,D$sexFemale)}
  if ('age' %in% covariates) {controls=cbind(controls,D$age)}
  if ('income' %in% covariates) {controls=cbind(controls,D$`income2: $25,000 to $50,000`,D$`income3: $50,001 to $75,000`,D$`income4: $75,001 to $100,000`,D$`income5: More than $100,001`)}
  if ('education' %in% covariates) {controls=cbind(controls,D$`educationSome College or Vocational`,D$`educationB.A. or B.S.`,D$`educationPost-graduate or Higher`)}
  if ('ethnicity' %in% covariates) {controls=cbind(controls,D$ethnicityBlack,D$`ethnicityNative American`,D$`ethnicityOther/Decline to State`,D$ethnicityWhite)}
  
  # new data to simulate over
  d.position.opp=cbind(controls,0,1,0,0,0,0,0,0)
  d.order.opp=cbind(controls,0,1,0,1,0,0,0,1)
  d.congress.opp=cbind(controls,0,1,1,0,0,1,0,0)
  
  d.position.cop=cbind(controls,1,0,0,0,0,0,0,0)
  d.order.cop=cbind(controls,1,0,0,1,0,0,1,0)
  d.congress.cop=cbind(controls,1,0,1,0,1,0,0,0)
  
  d.position.ind=cbind(controls,0,0,0,0,0,0,0,0)
  d.order.ind=cbind(controls,0,0,0,1,0,0,0,0)
  d.congress.ind=cbind(controls,0,0,1,0,0,0,0,0)
  
  # simulations
  p.pos.opp=apply((1/(1+exp(-as.matrix(d.position.opp)%*% t(betas)))),1,as.vector) 
  p_pos_opp=apply(p.pos.opp,1,mean)
  p.ord.opp=apply((1/(1+exp(-as.matrix(d.order.opp)%*% t(betas)))),1,as.vector) 
  p_ord_opp=apply(p.ord.opp,1,mean)
  p.cong.opp=apply((1/(1+exp(-as.matrix(d.congress.opp)%*% t(betas)))),1,as.vector) 
  p_cong_opp=apply(p.cong.opp,1,mean)
  
  p.pos.cop=apply((1/(1+exp(-as.matrix(d.position.cop)%*% t(betas)))),1,as.vector) 
  p_pos_cop=apply(p.pos.cop,1,mean)
  p.ord.cop=apply((1/(1+exp(-as.matrix(d.order.cop)%*% t(betas)))),1,as.vector) 
  p_ord_cop=apply(p.ord.cop,1,mean)
  p.cong.cop=apply((1/(1+exp(-as.matrix(d.congress.cop)%*% t(betas)))),1,as.vector) 
  p_cong_cop=apply(p.cong.cop,1,mean)
  
  p.pos.ind=apply((1/(1+exp(-as.matrix(d.position.ind)%*% t(betas)))),1,as.vector) 
  p_pos_ind=apply(p.pos.ind,1,mean)
  p.ord.ind=apply((1/(1+exp(-as.matrix(d.order.ind)%*% t(betas)))),1,as.vector) 
  p_ord_ind=apply(p.ord.ind,1,mean)
  p.cong.ind=apply((1/(1+exp(-as.matrix(d.congress.ind)%*% t(betas)))),1,as.vector) 
  p_cong_ind=apply(p.cong.ind,1,mean)
  
  # store
  results[1,1]<- quantile(p_ord_opp-p_pos_opp,.0083)
  results[1,2]<- quantile(p_ord_opp-p_pos_opp,.025)
  results[1,3]<- mean(p_ord_opp-p_pos_opp)
  results[1,4]<- quantile(p_ord_opp-p_pos_opp,.975)
  results[1,5]<- quantile(p_ord_opp-p_pos_opp,.9916)
  results[1,6]<- 'Executive Order'
  results[1,7]<- covariates[1]
  results[1,8]<- 'Opposition'
  results[2,1]<- quantile(p_cong_opp-p_pos_opp,.0083)
  results[2,2]<- quantile(p_cong_opp-p_pos_opp,.025)
  results[2,3]<- mean(p_cong_opp-p_pos_opp)
  results[2,4]<- quantile(p_cong_opp-p_pos_opp,.975)
  results[2,5]<- quantile(p_cong_opp-p_pos_opp,.9916) 
  results[2,6]<- 'Congress'
  results[2,7]<- covariates[1]
  results[2,8]<- 'Opposition'
  results[3,1]<- quantile(p_ord_opp-p_cong_opp,.0083)
  results[3,2]<- quantile(p_ord_opp-p_cong_opp,.025)
  results[3,3]<- mean(p_ord_opp-p_cong_opp)
  results[3,4]<- quantile(p_ord_opp-p_cong_opp,.975)
  results[3,5]<- quantile(p_ord_opp-p_cong_opp,.9916)
  results[3,6]<- 'Executive Order - Congress'
  results[3,7]<- covariates[1]
  results[3,8]<- 'Opposition'
  
  results[4,1]<- quantile(p_ord_cop-p_pos_cop,.0083)
  results[4,2]<- quantile(p_ord_cop-p_pos_cop,.025)
  results[4,3]<- mean(p_ord_cop-p_pos_cop)
  results[4,4]<- quantile(p_ord_cop-p_pos_cop,.975)
  results[4,5]<- quantile(p_ord_cop-p_pos_cop,.9916)
  results[4,6]<- 'Executive Order'
  results[4,7]<- covariates[1]
  results[4,8]<- 'Copartisans'
  results[5,1]<- quantile(p_cong_cop-p_pos_cop,.0083)
  results[5,2]<- quantile(p_cong_cop-p_pos_cop,.025)
  results[5,3]<- mean(p_cong_cop-p_pos_cop)
  results[5,4]<- quantile(p_cong_cop-p_pos_cop,.975)
  results[5,5]<- quantile(p_cong_cop-p_pos_cop,.9916) 
  results[5,6]<- 'Congress'
  results[5,7]<- covariates[1]
  results[5,8]<- 'Copartisans'
  results[6,1]<- quantile(p_ord_cop-p_cong_cop,.0083)
  results[6,2]<- quantile(p_ord_cop-p_cong_cop,.025)
  results[6,3]<- mean(p_ord_cop-p_cong_cop)
  results[6,4]<- quantile(p_ord_cop-p_cong_cop,.975)
  results[6,5]<- quantile(p_ord_cop-p_cong_cop,.9916)
  results[6,6]<- 'Executive Order - Congress'
  results[6,7]<- covariates[1]
  results[6,8]<- 'Copartisans'
  
  results[7,1]<- quantile(p_ord_ind-p_pos_ind,.0083)
  results[7,2]<- quantile(p_ord_ind-p_pos_ind,.025)
  results[7,3]<- mean(p_ord_ind-p_pos_ind)
  results[7,4]<- quantile(p_ord_ind-p_pos_ind,.975)
  results[7,5]<- quantile(p_ord_ind-p_pos_ind,.9916)
  results[7,6]<- 'Executive Order'
  results[7,7]<- covariates[1]
  results[7,8]<- 'Independents'
  results[8,1]<- quantile(p_cong_ind-p_pos_ind,.0083)
  results[8,2]<- quantile(p_cong_ind-p_pos_ind,.025)
  results[8,3]<- mean(p_cong_ind-p_pos_ind)
  results[8,4]<- quantile(p_cong_ind-p_pos_ind,.975)
  results[8,5]<- quantile(p_cong_ind-p_pos_ind,.9916) 
  results[8,6]<- 'Congress'
  results[8,7]<- covariates[1]
  results[8,8]<- 'Independents'
  results[9,1]<- quantile(p_ord_ind-p_cong_ind,.0083)
  results[9,2]<- quantile(p_ord_ind-p_cong_ind,.025)
  results[9,3]<- mean(p_ord_ind-p_cong_ind)
  results[9,4]<- quantile(p_ord_ind-p_cong_ind,.975)
  results[9,5]<- quantile(p_ord_ind-p_cong_ind,.9916)
  results[9,6]<- 'Executive Order - Congress'
  results[9,7]<- covariates[1]
  results[9,8]<- 'Independents'
  
  return(results)
}

outcomes=data.frame(lb.MC=numeric(54),lb.95=numeric(54),me=numeric(54),ub.95=numeric(54),ub.MC=numeric(54),comparison=character(54),dv=character(54),party=character(54))

outcomes[1:9,1:8]=public_comparisons_party(t3a.bin,dt)
outcomes[10:18,1:8]=public_comparisons_party(p3a.bin,dt)
outcomes[19:27,1:8]=public_comparisons_party(v3a,dt)
outcomes[28:36,1:8]=public_comparisons_party(c3a.bin,dt)
outcomes[37:45,1:8]=public_comparisons_party(g3a.bin,dt)
outcomes[46:54,1:8]=public_comparisons_party(l3a.bin,dt)

outcomes_overall$party='All Respondents'
final_outcomes=rbind(outcomes_overall,outcomes)

final_outcomes$lb.MC=as.numeric(final_outcomes$lb.MC)
final_outcomes$lb.95=as.numeric(final_outcomes$lb.95)
final_outcomes$me=as.numeric(final_outcomes$me)
final_outcomes$ub.95=as.numeric(final_outcomes$ub.95)
final_outcomes$ub.MC=as.numeric(final_outcomes$ub.MC)

final_outcomes$dv=as.factor(final_outcomes$dv)
levels(final_outcomes$dv)[levels(final_outcomes$dv)=='trumpvote.w1'] <- 'Trump\nVote'
levels(final_outcomes$dv)[levels(final_outcomes$dv)=='topic.approve.w1.bin'] <- 'Topic\nApproval'
levels(final_outcomes$dv)[levels(final_outcomes$dv)=='pres.approve.w1.bin'] <- 'Job\nApproval'
levels(final_outcomes$dv)[levels(final_outcomes$dv)=='manip.law.w1.bin'] <- 'Respects the\nLaw'
levels(final_outcomes$dv)[levels(final_outcomes$dv)=='manip.done.w1.bin'] <- 'Gets things\nDone'
levels(final_outcomes$dv)[levels(final_outcomes$dv)=='manip.cong.w1.bin'] <- 'Works with\nCongress'

other_outcomes=c('Works with\nCongress','Gets things\nDone','Respects the\nLaw')

# main h figure
f1=ggplot(final_outcomes[final_outcomes$comparison!='Executive Order - Congress' & !(final_outcomes$dv %in% other_outcomes),], aes(colour = comparison)) + 
  scale_colour_manual(name='...compared to position-taking.',values = c("#AA3377", "#004488")) +
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) + 
  geom_linerange(aes(x = factor(dv, 
                                levels = c("Trump\nVote", 
                                           "Job\nApproval", 
                                           "Topic\nApproval")),
                     ymin = lb.MC,
                     ymax = ub.MC),
                 lwd = 1.1, 
                 position = position_dodge(width = .5)) +
  geom_pointrange(aes(x = factor(dv, 
                                 levels = c("Trump\nVote", 
                                            "Job\nApproval", 
                                            "Topic\nApproval")),
                      y = me,
                      ymin = lb.95,
                      ymax = ub.95),
                  size=2.5,
                  lwd = 2,
                  position = position_dodge(width = .5),
                  shape = 21, 
                  fill = "WHITE") +
  geom_text(aes(x = factor(dv, levels = c("Trump\nVote", 
                                          "Job\nApproval", 
                                          "Topic\nApproval")),
                y = me,
                label=round(me,2)),
            size=3,
            position = position_dodge(width = .5)) +
  facet_wrap(~party,
             ncol=4,
             scales='free_x') +
  coord_flip() + 
  theme_bw() + 
  xlab('') + 
  ylab('') + 
  theme(text=element_text(family='Palatino'),
        legend.position='bottom',
        legend.title=element_text(size=9),
        strip.background = element_rect(
          color = "black", 
          fill = "#ffffff", 
          linewidth = 1.5, 
          linetype = "blank")) +
  guides(color = guide_legend(
    override.aes = list(shape = 19, size = 1),
    reverse = T,
    title.position = "right")) 
f1

ggsave("figures-tables/FIG-H1-1.pdf",
       f1,
       width = 7.2,
       height = 4,
       units = 'in')

# main h figure; other outcomes for appendix
f1_si = ggplot(final_outcomes[final_outcomes$comparison != 'Executive Order - Congress' &
                                final_outcomes$dv %in% other_outcomes, ], aes(colour = comparison)) +
  scale_colour_manual(name = '...compared to position-taking.', values = c("#AA3377", "#004488")) +
  geom_hline(yintercept = 0,
             colour = gray(1 / 2),
             lty = 2) +
  geom_linerange(aes(x = dv, ymin = lb.MC, ymax = ub.MC),
                 lwd = 1.1,
                 position = position_dodge(width = 1 / 2)) +
  geom_pointrange(
    aes(
      x = dv,
      y = me,
      ymin = lb.95,
      ymax = ub.95
    ),
    lwd = 2,
    size = 2.5,
    position = position_dodge(width = 1 / 2),
    shape = 21,
    fill = "WHITE"
  ) +
  geom_text(aes(x = dv, y = me, label = round(me, 2)),
            size = 3,
            position = position_dodge(width = 1 / 2)) +
  facet_wrap( ~ party, ncol = 4, scales = 'free_x') +
  coord_flip() + theme_bw() + xlab('') + ylab('') +
  theme(
    legend.position = 'bottom',
    legend.title = element_text(size = 9),
    strip.background = element_rect(
      color = "black",
      fill = "#ffffff",
      linewidth = 1.5,
      linetype = "blank"
    )
  ) +
  guides(color = guide_legend(
    override.aes = list(shape = 19, size = 1),
    reverse = T,
    title.position = "right"
  ))

ggsave("figures-tables/FIG-H1-A-1.pdf",
       f1_si,
       width = 7.2,
       height = 4,
       units = 'in')


# Save the tabular output used for the main text figure
stargazer(t3.bin,p3.bin,v3,c3.bin,g3.bin,l3.bin,digits=3,
          out = 'figures-tables/h1-tabular-output.txt',
          keep='condition',
          omit.stat = c('ll','aic'),
          omit.table.layout = 'n',
          font.size='footnotesize',
          column.sep.width='3pt',
          column.labels=c('\\shortstack{Topic\\\\Approval}', '\\shortstack{Job\\\\ Approval}', '\\shortstack{Trump\\\\ Vote}', '\\shortstack{Works w/\\\\ Congress}','\\shortstack{Gets things\\\\ Done}','\\shortstack{Cares about\\\\ Law}'),
          covariate.labels = c('Congress','Executive Order'),
          colnames = F, 
          dep.var.labels.include =F,
          model.numbers=F,
          title='{\\bf Mixed Evaluations of Unilateralism by the Public.} Reports logistic regression coefficients and conventional standard errors with binary dependent variables indicating either approval, voting for incumbent, or positive agreement with questions; simulated marginal effect estimates based on observed case approach from these models are reported in the leftmost panel of Figures \\ref{fig:h1} and \\ref{fig:h1a}; all models include topic and president factor variables, along with partisanship, age, income, sex, ethnicity, and education.',
          label='tab:h1',
          header=F)

stargazer(t3a.bin,p3a.bin,v3a,c3a.bin,g3a.bin,l3a.bin,digits=3,
          out = 'figures-tables/h1-tabular-output-2.txt',
          keep=c('condition','copartisanship'),
          omit.stat = c('ll','aic'),
          omit.table.layout = 'n',
          font.size='tiny',
          column.labels=c('\\shortstack{Topic\\\\Approval}', '\\shortstack{Job\\\\ Approval}', '\\shortstack{Trump\\\\ Vote}', '\\shortstack{Works w/\\\\ Congress}','\\shortstack{Gets things\\\\ Done}','\\shortstack{Cares about\\\\ Law}'),
          column.sep.width='3pt',
          colnames = F,
          covariate.labels=c('Copartisan','Opposition','Congress','Exec. Order','Copartisan X Congress','Opposition X Congress','Copartisan X Exec. Order','Opposition X Exec. Order'),
          dep.var.labels.include =F,
          model.numbers=F,
          title='{\\bf Mixed Evaluations of Unilateralism by Copartisans, the Opposition, and Independents.} Reports logistic regression coefficients and conventional standard errors with binary dependent variables indicating either approval, voting for incumbent, or positive agreement with questions; simulated marginal effect estimates based on observed case approach from these models are reported in the right three panels of Figures \\ref{fig:h1} and \\ref{fig:h1a}; all models include topic and president factor variables, along with age, income, sex, ethnicity, and education.',
          label='tab:h1a',
          header=F)


stargazer(t3, p3, w3, g3, l3, 
          digits = 3, 
          out = 'figures-tables/h1-tabular-output-ordered.txt',
          column.labels=c('\\shortstack{Topic\\\\Approval}', '\\shortstack{Job\\\\ Approval}', '\\shortstack{Works w/\\\\ Congress}','\\shortstack{Gets things\\\\ Done}','\\shortstack{Cares about\\\\ Law}'),
          column.sep.width='3pt',
          covariate.labels = c("Congress", "Executive Order"),
          omit.stat = c('ll', 'aic'),
          omit.table.layout = 'n',
          font.size = 'tiny',
          keep=c('condition'),
          colnames = F,
          header = F,
          dep.var.labels.include =F,
          model.numbers=F,
          title='{\\bf Mixed Evaluations of Unilateralism by the Public (Ordered Logits).} Reports ordered logistic regression coefficients and conventional standard errors with 7-point Likert dependent variables indicating ``Strongly disagree,`` ``Disagree,`` ``Somewhat disagree,`` ``Neither agree nor disagree,`` ``Somewhat agree,`` ``Agree,`` ``Strongly agree;`` all models include topic and president factor variables, along with partisanship, age, income, sex, ethnicity, and education.',
          label='tab:h1_ordered')


stargazer(t3a, p3a, w3a, g3a, l3a, 
          digits = 3, 
          out = 'figures-tables/h1-tabular-output-ordered-2.txt',
          column.labels=c('\\shortstack{Topic\\\\Approval}', '\\shortstack{Job\\\\ Approval}', '\\shortstack{Works w/\\\\ Congress}','\\shortstack{Gets things\\\\ Done}','\\shortstack{Cares about\\\\ Law}'),
          column.sep.width='3pt',
          omit.stat = c('ll', 'aic'),
          omit.table.layout = 'n',
          font.size = 'tiny',
          covariate.labels=c('Copartisan','Opposition','Congress','Exec. Order','Copartisan X Congress','Opposition X Congress','Copartisan X Exec. Order','Opposition X Exec. Order'),
          keep=c('condition', 'copartisanship'),
          colnames = F,
          header = F,
          dep.var.labels.include =F,
          model.numbers=F,
          title='{\\bf Mixed Evaluations of Unilateralism by Copartisans, the Opposition, and Independents (Ordered Logits).} Reports ordered logistic regression coefficients and conventional standard errors with 7-point Likert dependent variables indicating ``Strongly disagree,`` ``Disagree,`` ``Somewhat disagree,`` ``Neither agree nor disagree,`` ``Somewhat agree,`` ``Agree,`` ``Strongly agree;`` all models include topic and president factor variables, along with age, income, sex, ethnicity, and education.',
          label='tab:h1a_ordered')


################################################
## figure 1 w/same-party 2020 vote (R1 and R2) #
################################################

# SAME-PARTY VOTE
dt$sameparty.w1 <- 0
dt$sameparty.w1[dt$president.w1=='Obama' & dt$votechoice.w1=='Joseph R. Biden']=1
dt$sameparty.w1[dt$president.w1=='Trump' & dt$votechoice.w1=='Donald J. Trump']=1

# TOPIC APPROVAL
t1=clm(topic.approve.w1~topic.w1+president.w1+condition.w1,data=dt,family='logit')
t2=clm(topic.approve.w1~topic.w1+president.w1+partisanship+sex+age+income+ethnicity+condition.w1,data=dt,family='logit')
t3=clm(topic.approve.w1~topic.w1+president.w1+partisanship+sex+age+income+education+ethnicity+condition.w1,data=dt,family='logit')

# PRESIDENTIAL APPROVAL
p1=clm(pres.approve.w1~topic.w1+president.w1+condition.w1,data=dt,family='logit')
p2=clm(pres.approve.w1~topic.w1+president.w1+partisanship+sex+age+income+ethnicity+condition.w1,data=dt,family='logit')
p3=clm(pres.approve.w1~topic.w1+president.w1+partisanship+sex+age+income+education+ethnicity+condition.w1,data=dt,family='logit')

# TRUMP VOTE
v1=glm(sameparty.w1~topic.w1+president.w1+condition.w1,data=dt,family='binomial')
v2=glm(sameparty.w1~topic.w1+president.w1+partisanship+sex+age+income+ethnicity+condition.w1,data=dt,family='binomial')
v3=glm(sameparty.w1~topic.w1+president.w1+partisanship+sex+age+income+education+ethnicity+condition.w1,data=dt,family='binomial')

# Estimates for main figure

# re-run as logits
t3.bin=glm(topic.approve.w1.bin~topic.w1+president.w1+partisanship+sex+age+income+education+ethnicity+condition.w1,data=dt,family='binomial')
p3.bin=glm(pres.approve.w1.bin~topic.w1+president.w1+partisanship+sex+age+income+education+ethnicity+condition.w1,data=dt,family='binomial')

outcomes=matrix(data=NA, ncol=7,nrow=9)

public_comparisons=function(model,data,simulations=1000,seed=6658) {
  set.seed(seed)
  betas=slot(sim(model,simulations),'coef')
  covariates=names(model$model)
  results=matrix(data=NA, ncol=7,nrow=3)
  
  # add columns for factor coefs
  D=data[complete.cases(data[,covariates]),]
  D=cbind(D,model.matrix(~topic.w1-1,data=D))
  D=cbind(D,model.matrix(~president.w1-1,data=D))
  if ('partisanship' %in% covariates) {D=cbind(D,model.matrix(~partisanship-1,data=D))}
  if ('sex' %in% covariates) {D=cbind(D,model.matrix(~sex-1,data=D))}
  if ('income' %in% covariates) {D=cbind(D,model.matrix(~income-1,data=D))}
  if ('education' %in% covariates) {D=cbind(D,model.matrix(~education-1,data=D))}
  if ('ethnicity' %in% covariates) {D=cbind(D,model.matrix(~ethnicity-1,data=D))}
  D=cbind(D,model.matrix(~condition.w1-1,data=D))
  
  # intercept + covariates 
  controls=cbind(1,D$`topic.w1contractor pay`,D$`topic.w1farm payments`,D$`topic.w1foreign aid and abortion`,D$`topic.w1foreign workers`,D$`topic.w1gun violence research`,D$topic.w1LGBT,D$`topic.w1policing and military surplus`,D$`topic.w1public lands`,D$`topic.w1Russian sanctions`,D$`topic.w1student loans`,D$topic.w1trade,D$`topic.w1water rules`,D$`topic.w1weapon sales to Saudi Arabia`,D$topic.w1wildlife,D$president.w1Trump)
  if ('partisanship' %in% covariates) {controls=cbind(controls,D$partisanshipDemocrat,D$partisanshipRepublican)}
  if ('sex' %in% covariates) {controls=cbind(controls,D$sexFemale)}
  if ('age' %in% covariates) {controls=cbind(controls,D$age)}
  if ('income' %in% covariates) {controls=cbind(controls,D$`income2: $25,000 to $50,000`,D$`income3: $50,001 to $75,000`,D$`income4: $75,001 to $100,000`,D$`income5: More than $100,001`)}
  if ('education' %in% covariates) {controls=cbind(controls,D$`educationSome College or Vocational`,D$`educationB.A. or B.S.`,D$`educationPost-graduate or Higher`)}
  if ('ethnicity' %in% covariates) {controls=cbind(controls,D$ethnicityBlack,D$`ethnicityNative American`,D$`ethnicityOther/Decline to State`,D$ethnicityWhite)}
  
  # new data to simulate over
  d.position=cbind(controls,0,0)
  d.order=cbind(controls,0,1)
  d.congress=cbind(controls,1,0)
  
  # simulations
  p.pos=apply((1/(1+exp(-as.matrix(d.position)%*% t(betas)))),1,as.vector) 
  p_pos=apply(p.pos,1,mean)
  p.ord=apply((1/(1+exp(-as.matrix(d.order)%*% t(betas)))),1,as.vector) 
  p_ord=apply(p.ord,1,mean)
  p.cong=apply((1/(1+exp(-as.matrix(d.congress)%*% t(betas)))),1,as.vector) 
  p_cong=apply(p.cong,1,mean)
  
  # store
  results[1,1]<- quantile(p_ord-p_pos,.0083)
  results[1,2]<- quantile(p_ord-p_pos,.025)
  results[1,3]<- mean(p_ord-p_pos)
  results[1,4]<- quantile(p_ord-p_pos,.975)
  results[1,5]<- quantile(p_ord-p_pos,.9916)
  results[1,6]<- 'Executive Order'
  results[1,7]<- covariates[1]
  results[2,1]<- quantile(p_cong-p_pos,.0083)
  results[2,2]<- quantile(p_cong-p_pos,.025)
  results[2,3]<- mean(p_cong-p_pos)
  results[2,4]<- quantile(p_cong-p_pos,.975)
  results[2,5]<- quantile(p_cong-p_pos,.9916) 
  results[2,6]<- 'Congress'
  results[2,7]<- covariates[1]
  results[3,1]<- quantile(p_ord-p_cong,.0083)
  results[3,2]<- quantile(p_ord-p_cong,.025)
  results[3,3]<- mean(p_ord-p_cong)
  results[3,4]<- quantile(p_ord-p_cong,.975)
  results[3,5]<- quantile(p_ord-p_cong,.9916)
  results[3,6]<- 'Executive Order - Congress'
  results[3,7]<- covariates[1]
  
  return(results)
}

outcomes=data.frame(lb.MC=numeric(9),lb.95=numeric(9),me=numeric(9),ub.95=numeric(9),ub.MC=numeric(9),comparison=character(9),dv=character(9))
outcomes[1:3,1:7]=public_comparisons(t3.bin,dt)
outcomes[4:6,1:7]=public_comparisons(p3.bin,dt)
outcomes[7:9,1:7]=public_comparisons(v3,dt)
outcomes_overall=outcomes

# hypothesis 1A: break down co-partisanship

# TOPIC APPROVAL
t1a=clm(topic.approve.w1~topic.w1+president.w1+copartisanship*condition.w1,data=dt,family='logit')
t2a=clm(topic.approve.w1~topic.w1+president.w1+sex+age+income+ethnicity+copartisanship*condition.w1,data=dt,family='logit')
t3a=clm(topic.approve.w1~topic.w1+president.w1+sex+age+income+education+ethnicity+copartisanship*condition.w1,data=dt,family='logit')

# PRESIDENTIAL APPROVAL
p1a=clm(pres.approve.w1~topic.w1+president.w1+copartisanship*condition.w1,data=dt,family='logit')
p2a=clm(pres.approve.w1~topic.w1+president.w1+sex+age+income+ethnicity+copartisanship*condition.w1,data=dt,family='logit')
p3a=clm(pres.approve.w1~topic.w1+president.w1+sex+age+income+education+ethnicity+copartisanship*condition.w1,data=dt,family='logit')

# SAME-PARTY VOTE
v1a=glm(sameparty.w1~topic.w1+president.w1+copartisanship*condition.w1,data=dt,family='binomial')
v2a=glm(sameparty.w1~topic.w1+president.w1+sex+age+income+ethnicity+copartisanship*condition.w1,data=dt,family='binomial')
v3a=glm(sameparty.w1~topic.w1+president.w1+sex+age+income+education+ethnicity+copartisanship*condition.w1,data=dt,family='binomial')

# re-run as logits
t3a.bin=glm(topic.approve.w1.bin~topic.w1+president.w1+sex+age+income+education+ethnicity+copartisanship*condition.w1,data=dt,family='binomial')
p3a.bin=glm(pres.approve.w1.bin~topic.w1+president.w1+sex+age+income+education+ethnicity+copartisanship*condition.w1,data=dt,family='binomial')
v3a=glm(sameparty.w1~topic.w1+president.w1+sex+age+income+education+ethnicity+copartisanship*condition.w1,data=dt,family='binomial')

public_comparisons_party=function(model,data,simulations=1000,seed=6658) {
  set.seed(seed)
  betas=slot(sim(model,simulations),'coef')
  covariates=names(model$model)
  results=matrix(data=NA, ncol=8,nrow=9)
  
  # add columns for factor coefs
  D=data[complete.cases(data[,covariates]),]
  D=cbind(D,model.matrix(~topic.w1-1,data=D))
  D=cbind(D,model.matrix(~president.w1-1,data=D))
  if ('copartisanship' %in% covariates) {D=cbind(D,model.matrix(~copartisanship-1,data=D))}
  if ('sex' %in% covariates) {D=cbind(D,model.matrix(~sex-1,data=D))}
  if ('income' %in% covariates) {D=cbind(D,model.matrix(~income-1,data=D))}
  if ('education' %in% covariates) {D=cbind(D,model.matrix(~education-1,data=D))}
  if ('ethnicity' %in% covariates) {D=cbind(D,model.matrix(~ethnicity-1,data=D))}
  D=cbind(D,model.matrix(~condition.w1-1,data=D))
  
  # intercept + covariates 
  controls=cbind(1,D$`topic.w1contractor pay`,D$`topic.w1farm payments`,D$`topic.w1foreign aid and abortion`,D$`topic.w1foreign workers`,D$`topic.w1gun violence research`,D$topic.w1LGBT,D$`topic.w1policing and military surplus`,D$`topic.w1public lands`,D$`topic.w1Russian sanctions`,D$`topic.w1student loans`,D$topic.w1trade,D$`topic.w1water rules`,D$`topic.w1weapon sales to Saudi Arabia`,D$topic.w1wildlife,D$president.w1Trump)
  if ('sex' %in% covariates) {controls=cbind(controls,D$sexFemale)}
  if ('age' %in% covariates) {controls=cbind(controls,D$age)}
  if ('income' %in% covariates) {controls=cbind(controls,D$`income2: $25,000 to $50,000`,D$`income3: $50,001 to $75,000`,D$`income4: $75,001 to $100,000`,D$`income5: More than $100,001`)}
  if ('education' %in% covariates) {controls=cbind(controls,D$`educationSome College or Vocational`,D$`educationB.A. or B.S.`,D$`educationPost-graduate or Higher`)}
  if ('ethnicity' %in% covariates) {controls=cbind(controls,D$ethnicityBlack,D$`ethnicityNative American`,D$`ethnicityOther/Decline to State`,D$ethnicityWhite)}
  
  # new data to simulate over
  d.position.opp=cbind(controls,0,1,0,0,0,0,0,0)
  d.order.opp=cbind(controls,0,1,0,1,0,0,0,1)
  d.congress.opp=cbind(controls,0,1,1,0,0,1,0,0)
  
  d.position.cop=cbind(controls,1,0,0,0,0,0,0,0)
  d.order.cop=cbind(controls,1,0,0,1,0,0,1,0)
  d.congress.cop=cbind(controls,1,0,1,0,1,0,0,0)
  
  d.position.ind=cbind(controls,0,0,0,0,0,0,0,0)
  d.order.ind=cbind(controls,0,0,0,1,0,0,0,0)
  d.congress.ind=cbind(controls,0,0,1,0,0,0,0,0)
  
  # simulations
  p.pos.opp=apply((1/(1+exp(-as.matrix(d.position.opp)%*% t(betas)))),1,as.vector) 
  p_pos_opp=apply(p.pos.opp,1,mean)
  p.ord.opp=apply((1/(1+exp(-as.matrix(d.order.opp)%*% t(betas)))),1,as.vector) 
  p_ord_opp=apply(p.ord.opp,1,mean)
  p.cong.opp=apply((1/(1+exp(-as.matrix(d.congress.opp)%*% t(betas)))),1,as.vector) 
  p_cong_opp=apply(p.cong.opp,1,mean)
  
  p.pos.cop=apply((1/(1+exp(-as.matrix(d.position.cop)%*% t(betas)))),1,as.vector) 
  p_pos_cop=apply(p.pos.cop,1,mean)
  p.ord.cop=apply((1/(1+exp(-as.matrix(d.order.cop)%*% t(betas)))),1,as.vector) 
  p_ord_cop=apply(p.ord.cop,1,mean)
  p.cong.cop=apply((1/(1+exp(-as.matrix(d.congress.cop)%*% t(betas)))),1,as.vector) 
  p_cong_cop=apply(p.cong.cop,1,mean)
  
  p.pos.ind=apply((1/(1+exp(-as.matrix(d.position.ind)%*% t(betas)))),1,as.vector) 
  p_pos_ind=apply(p.pos.ind,1,mean)
  p.ord.ind=apply((1/(1+exp(-as.matrix(d.order.ind)%*% t(betas)))),1,as.vector) 
  p_ord_ind=apply(p.ord.ind,1,mean)
  p.cong.ind=apply((1/(1+exp(-as.matrix(d.congress.ind)%*% t(betas)))),1,as.vector) 
  p_cong_ind=apply(p.cong.ind,1,mean)
  
  # store
  results[1,1]<- quantile(p_ord_opp-p_pos_opp,.0083)
  results[1,2]<- quantile(p_ord_opp-p_pos_opp,.025)
  results[1,3]<- mean(p_ord_opp-p_pos_opp)
  results[1,4]<- quantile(p_ord_opp-p_pos_opp,.975)
  results[1,5]<- quantile(p_ord_opp-p_pos_opp,.9916)
  results[1,6]<- 'Executive Order'
  results[1,7]<- covariates[1]
  results[1,8]<- 'Opposition'
  results[2,1]<- quantile(p_cong_opp-p_pos_opp,.0083)
  results[2,2]<- quantile(p_cong_opp-p_pos_opp,.025)
  results[2,3]<- mean(p_cong_opp-p_pos_opp)
  results[2,4]<- quantile(p_cong_opp-p_pos_opp,.975)
  results[2,5]<- quantile(p_cong_opp-p_pos_opp,.9916) 
  results[2,6]<- 'Congress'
  results[2,7]<- covariates[1]
  results[2,8]<- 'Opposition'
  results[3,1]<- quantile(p_ord_opp-p_cong_opp,.0083)
  results[3,2]<- quantile(p_ord_opp-p_cong_opp,.025)
  results[3,3]<- mean(p_ord_opp-p_cong_opp)
  results[3,4]<- quantile(p_ord_opp-p_cong_opp,.975)
  results[3,5]<- quantile(p_ord_opp-p_cong_opp,.9916)
  results[3,6]<- 'Executive Order - Congress'
  results[3,7]<- covariates[1]
  results[3,8]<- 'Opposition'
  
  results[4,1]<- quantile(p_ord_cop-p_pos_cop,.0083)
  results[4,2]<- quantile(p_ord_cop-p_pos_cop,.025)
  results[4,3]<- mean(p_ord_cop-p_pos_cop)
  results[4,4]<- quantile(p_ord_cop-p_pos_cop,.975)
  results[4,5]<- quantile(p_ord_cop-p_pos_cop,.9916)
  results[4,6]<- 'Executive Order'
  results[4,7]<- covariates[1]
  results[4,8]<- 'Copartisans'
  results[5,1]<- quantile(p_cong_cop-p_pos_cop,.0083)
  results[5,2]<- quantile(p_cong_cop-p_pos_cop,.025)
  results[5,3]<- mean(p_cong_cop-p_pos_cop)
  results[5,4]<- quantile(p_cong_cop-p_pos_cop,.975)
  results[5,5]<- quantile(p_cong_cop-p_pos_cop,.9916) 
  results[5,6]<- 'Congress'
  results[5,7]<- covariates[1]
  results[5,8]<- 'Copartisans'
  results[6,1]<- quantile(p_ord_cop-p_cong_cop,.0083)
  results[6,2]<- quantile(p_ord_cop-p_cong_cop,.025)
  results[6,3]<- mean(p_ord_cop-p_cong_cop)
  results[6,4]<- quantile(p_ord_cop-p_cong_cop,.975)
  results[6,5]<- quantile(p_ord_cop-p_cong_cop,.9916)
  results[6,6]<- 'Executive Order - Congress'
  results[6,7]<- covariates[1]
  results[6,8]<- 'Copartisans'
  
  results[7,1]<- quantile(p_ord_ind-p_pos_ind,.0083)
  results[7,2]<- quantile(p_ord_ind-p_pos_ind,.025)
  results[7,3]<- mean(p_ord_ind-p_pos_ind)
  results[7,4]<- quantile(p_ord_ind-p_pos_ind,.975)
  results[7,5]<- quantile(p_ord_ind-p_pos_ind,.9916)
  results[7,6]<- 'Executive Order'
  results[7,7]<- covariates[1]
  results[7,8]<- 'Independents'
  results[8,1]<- quantile(p_cong_ind-p_pos_ind,.0083)
  results[8,2]<- quantile(p_cong_ind-p_pos_ind,.025)
  results[8,3]<- mean(p_cong_ind-p_pos_ind)
  results[8,4]<- quantile(p_cong_ind-p_pos_ind,.975)
  results[8,5]<- quantile(p_cong_ind-p_pos_ind,.9916) 
  results[8,6]<- 'Congress'
  results[8,7]<- covariates[1]
  results[8,8]<- 'Independents'
  results[9,1]<- quantile(p_ord_ind-p_cong_ind,.0083)
  results[9,2]<- quantile(p_ord_ind-p_cong_ind,.025)
  results[9,3]<- mean(p_ord_ind-p_cong_ind)
  results[9,4]<- quantile(p_ord_ind-p_cong_ind,.975)
  results[9,5]<- quantile(p_ord_ind-p_cong_ind,.9916)
  results[9,6]<- 'Executive Order - Congress'
  results[9,7]<- covariates[1]
  results[9,8]<- 'Independents'
  
  return(results)
}

outcomes=data.frame(lb.MC=numeric(27),lb.95=numeric(27),me=numeric(27),ub.95=numeric(27),ub.MC=numeric(27),comparison=character(27),dv=character(27),party=character(27))

outcomes[1:9,1:8]=public_comparisons_party(t3a.bin,dt)
outcomes[10:18,1:8]=public_comparisons_party(p3a.bin,dt)
outcomes[19:27,1:8]=public_comparisons_party(v3a,dt)

outcomes_overall$party='All Respondents'
final_outcomes=rbind(outcomes_overall,outcomes)

final_outcomes$lb.MC=as.numeric(final_outcomes$lb.MC)
final_outcomes$lb.95=as.numeric(final_outcomes$lb.95)
final_outcomes$me=as.numeric(final_outcomes$me)
final_outcomes$ub.95=as.numeric(final_outcomes$ub.95)
final_outcomes$ub.MC=as.numeric(final_outcomes$ub.MC)

final_outcomes$dv=as.factor(final_outcomes$dv)
levels(final_outcomes$dv)[levels(final_outcomes$dv)=='sameparty.w1'] <- 'Same-Party\nVote'
levels(final_outcomes$dv)[levels(final_outcomes$dv)=='topic.approve.w1.bin'] <- 'Topic\nApproval'
levels(final_outcomes$dv)[levels(final_outcomes$dv)=='pres.approve.w1.bin'] <- 'Job\nApproval'

# main h figure
f1_b = ggplot(final_outcomes[final_outcomes$comparison != 'Executive Order - Congress' &
                               !(final_outcomes$dv %in% other_outcomes), ], aes(colour = comparison)) +
  scale_colour_manual(name = '...compared to position-taking.', values = c("#AA3377", "#004488")) +
  geom_hline(yintercept = 0,
             colour = gray(1 / 2),
             lty = 2) +
  geom_linerange(
    aes(
      x = factor(
        dv,
        levels = c('Same-Party\nVote', 'Job\nApproval', 'Topic\nApproval')
      ),
      ymin = lb.MC,
      ymax = ub.MC
    ),
    lwd = 1.1,
    position = position_dodge(width = .5)
  ) +
  geom_pointrange(
    aes(
      x = factor(
        dv,
        levels = c('Same-Party\nVote', 'Job\nApproval', 'Topic\nApproval')
      ),
      y = me,
      ymin = lb.95,
      ymax = ub.95
    ),
    size = 2.5,
    lwd = 2,
    position = position_dodge(width = .5),
    shape = 21,
    fill = "WHITE"
  ) +
  geom_text(aes(
    x = factor(
      dv,
      levels = c('Same-Party\nVote', 'Job\nApproval', 'Topic\nApproval')
    ),
    y = me,
    label = round(me, 2)
  ),
  size = 3,
  position = position_dodge(width = .5)) +
  facet_wrap( ~ party, ncol = 4, scales = 'free_x') +
  coord_flip() + theme_bw() + xlab('') + ylab('') +
  theme(
    text = element_text(family = 'Palatino'),
    legend.position = 'bottom',
    legend.title = element_text(size = 9),
    strip.background = element_rect(
      color = "black",
      fill = "#ffffff",
      linewidth = 1.5,
      linetype = "blank"
    )
  ) +
  guides(color = guide_legend(
    override.aes = list(shape = 19, size = 1),
    reverse = T,
    title.position = "right"
  ))

ggsave("figures-tables/sameparty_f2.pdf",
       f1_b,
       width = 7.2,
       height = 4,
       units = 'in')


################################################
## figure 3 w/same-party 2020 vote (R1 and R2) #
################################################

dt$sameparty.w2 <- 0
dt$sameparty.w2[is.na(dt$president.w2)] <- NA
dt$sameparty.w2[dt$president.w2=='Obama' & dt$votechoice.w2=='Joseph R. Biden']=1
dt$sameparty.w2[dt$president.w2=='Trump' & dt$votechoice.w2=='Donald J. Trump']=1


# TOPIC APPROVAL
t1_c=clm(topic.approve.w2~topic.w2+president.w2+condition.w2+topic.approve.w1+outcome,data=dt,family='logit')
t2_c=clm(topic.approve.w2~topic.w2+president.w2+condition.w2+topic.approve.w1+partisanship+sex+age+income+ethnicity+outcome,data=dt,family='logit')
t3_c=clm(topic.approve.w2~topic.w2+president.w2+condition.w2+topic.approve.w1+partisanship+sex+age+income+ethnicity+education+outcome,data=dt,family='logit')

# PRESIDENTIAL APPROVAL
p1_c=clm(pres.approve.w2~topic.w2+president.w2+condition.w2+pres.approve.w1+outcome,data=dt,family='logit')
p2_c=clm(pres.approve.w2~topic.w2+president.w2+condition.w2+pres.approve.w1+partisanship+sex+age+income+ethnicity+outcome,data=dt,family='logit')
p3_c=clm(pres.approve.w2~topic.w2+president.w2+condition.w2+pres.approve.w1+partisanship+sex+age+income+ethnicity+education+outcome,data=dt,family='logit')

# TRUMP VOTE
v1_c=glm(sameparty.w2~topic.w2+president.w2+condition.w2+sameparty.w1+outcome,data=dt,family='binomial')
v2_c=glm(sameparty.w2~topic.w2+president.w2+condition.w2+sameparty.w1+partisanship+sex+age+income+ethnicity+outcome,data=dt,family='binomial')
v3_c=glm(sameparty.w2~topic.w2+president.w2+condition.w2+sameparty.w1+partisanship+sex+age+income+ethnicity+education+outcome,data=dt,family='binomial')

# rerun models with binary DV
t3_c.bin=glm(topic.approve.w2.bin~topic.w2+president.w2+condition.w2+topic.approve.w1+partisanship+sex+age+income+ethnicity+education+outcome,data=dt,family='binomial')
p3_c.bin=glm(pres.approve.w2.bin~topic.w2+president.w2+condition.w2+pres.approve.w1+partisanship+sex+age+income+ethnicity+education+outcome,data=dt,family='binomial')

# function for simulating marginal changes
outcome_comparisons=function(model,data,simulations=1000,seed=6658) {
  set.seed(seed)
  betas=slot(sim(model,simulations),'coef')
  covariates=names(model$model)
  results=matrix(data=NA, ncol=7,nrow=1)
  
  # add columns for factor coefs
  D=data[complete.cases(data[,covariates]),]
  D=cbind(D,model.matrix(~topic.w2-1,data=D))
  D=cbind(D,model.matrix(~president.w2-1,data=D))
  if ('topic.approve.w1' %in% covariates) {D=cbind(D,model.matrix(~topic.approve.w1-1,data=D))}
  if ('pres.approve.w1' %in% covariates) {D=cbind(D,model.matrix(~pres.approve.w1-1,data=D))}
  if ('manip.done.w1' %in% covariates) {D=cbind(D,model.matrix(~manip.done.w1-1,data=D))}
  if ('partisanship' %in% covariates) {D=cbind(D,model.matrix(~partisanship-1,data=D))}
  if ('sex' %in% covariates) {D=cbind(D,model.matrix(~sex-1,data=D))}
  if ('income' %in% covariates) {D=cbind(D,model.matrix(~income-1,data=D))}
  if ('ethnicity' %in% covariates) {D=cbind(D,model.matrix(~ethnicity-1,data=D))}
  if ('education' %in% covariates) {D=cbind(D,model.matrix(~education-1,data=D))}
  D=cbind(D,model.matrix(~condition.w2-1,data=D))
  
  # intercept + covariates 
  controls=cbind(1,D$`topic.w2contractor pay`,D$`topic.w2farm payments`,D$`topic.w2foreign aid and abortion`,D$`topic.w2foreign workers`,D$`topic.w2gun violence research`,D$topic.w2LGBT,D$`topic.w2policing and military surplus`,D$`topic.w2public lands`,D$`topic.w2Russian sanctions`,D$`topic.w2student loans`,D$topic.w2trade,D$`topic.w2water rules`,D$`topic.w2weapon sales to Saudi Arabia`,D$topic.w2wildlife,D$president.w2Trump,D$condition.w2congress,D$condition.w2order)
  if ('topic.approve.w1' %in% covariates) {controls=cbind(controls,D$topic.approve.w1Disapprove,D$`topic.approve.w1Somewhat disapprove`,D$`topic.approve.w1Neither approve or disapprove`,D$`topic.approve.w1Somewhat approve`,D$`topic.approve.w1Approve`,D$`topic.approve.w1Strongly approve`)}
  if ('pres.approve.w1' %in% covariates) {controls=cbind(controls,D$pres.approve.w1Disapprove,D$`pres.approve.w1Somewhat disapprove`,D$`pres.approve.w1Neither approve or disapprove`,D$`pres.approve.w1Somewhat approve`,D$`pres.approve.w1Approve`,D$`pres.approve.w1Strongly approve`)}
  if ('manip.done.w1' %in% covariates) {controls=cbind(controls,D$manip.done.w1Disagree,D$`manip.done.w1Somewhat disagree`,D$`manip.done.w1Neither agree nor disagree`,D$`manip.done.w1Somewhat agree`,D$`manip.done.w1Agree`,D$`manip.done.w1Strongly agree`)}
  if ('sameparty.w1' %in% covariates) {controls=cbind(controls,D$sameparty.w1)}
  if ('partisanship' %in% covariates) {controls=cbind(controls,D$partisanshipDemocrat,D$partisanshipRepublican)}
  if ('sex' %in% covariates) {controls=cbind(controls,D$sexFemale)}
  if ('age' %in% covariates) {controls=cbind(controls,D$age)}
  if ('income' %in% covariates) {controls=cbind(controls,D$`income2: $25,000 to $50,000`,D$`income3: $50,001 to $75,000`,D$`income4: $75,001 to $100,000`,D$`income5: More than $100,001`)}
  if ('ethnicity' %in% covariates) {controls=cbind(controls,D$ethnicityBlack,D$`ethnicityNative American`,D$`ethnicityOther/Decline to State`,D$ethnicityWhite)}
  if ('education' %in% covariates) {controls=cbind(controls,D$`educationSome College or Vocational`,D$`educationB.A. or B.S.`,D$`educationPost-graduate or Higher`)}
  
  # new data to simulate over
  d.success=cbind(controls,0)
  d.fail=cbind(controls,1)
  
  # simulations
  p.win=apply((1/(1+exp(-as.matrix(d.success)%*% t(betas)))),1,as.vector) 
  p_win=apply(p.win,1,mean)
  p.fail=apply((1/(1+exp(-as.matrix(d.fail)%*% t(betas)))),1,as.vector) 
  p_fail=apply(p.fail,1,mean)
  
  # store
  results[1,1]<- quantile(p_fail-p_win,.0083)
  results[1,2]<- quantile(p_fail-p_win,.025)
  results[1,3]<- mean(p_fail-p_win)
  results[1,4]<- quantile(p_fail-p_win,.975)
  results[1,5]<- quantile(p_fail-p_win,.9916)
  results[1,6]<- 'All Respondents'
  results[1,7]<- covariates[1]
  
  return(results)
}

outcomes=data.frame(lb.MC=numeric(3),lb.95=numeric(3),me=numeric(3),ub.95=numeric(3),ub.MC=numeric(3),grp=character(3),dv=character(3))
outcomes[1,1:7]=outcome_comparisons(t3_c.bin,dt)
outcomes[2,1:7]=outcome_comparisons(p3_c.bin,dt)
outcomes[3,1:7]=outcome_comparisons(v3_c,dt)

# copartisanship interactions

# TOPIC APPROVAL
t1a_c=clm(topic.approve.w2~topic.w2+president.w2+condition.w2+topic.approve.w1+copartisanship*outcome,data=dt,family='logit')
t2a_c=clm(topic.approve.w2~topic.w2+president.w2+condition.w2+topic.approve.w1+sex+age+income+ethnicity+copartisanship*outcome,data=dt,family='logit')
t3a_c=clm(topic.approve.w2~topic.w2+president.w2+condition.w2+topic.approve.w1+sex+age+income+ethnicity+education+copartisanship*outcome,data=dt,family='logit')

# JOB APPROVAL
p1a_c=clm(pres.approve.w2~topic.w2+president.w2+condition.w2+pres.approve.w1+copartisanship*outcome,data=dt,family='logit')
p2a_c=clm(pres.approve.w2~topic.w2+president.w2+condition.w2+pres.approve.w1+sex+age+income+ethnicity+copartisanship*outcome,data=dt,family='logit')
p3a_c=clm(pres.approve.w2~topic.w2+president.w2+condition.w2+pres.approve.w1+sex+age+income+ethnicity+education+copartisanship*outcome,data=dt,family='logit')

# TRUMP VOTE
v1a_c=glm(sameparty.w2~topic.w2+president.w2+condition.w2+sameparty.w1+copartisanship*outcome,data=dt,family='binomial')
v2a_c=glm(sameparty.w2~topic.w2+president.w2+condition.w2+sameparty.w1+sex+age+income+ethnicity+copartisanship*outcome,data=dt,family='binomial')
v3a_c=glm(sameparty.w2~topic.w2+president.w2+condition.w2+sameparty.w1+sex+age+income+ethnicity+education+copartisanship*outcome,data=dt,family='binomial')

# rerun with binary DV for simulations
t3a_c.bin=glm(topic.approve.w2.bin~topic.w2+president.w2+condition.w2+topic.approve.w1+sex+age+income+ethnicity+copartisanship*outcome,data=dt,family='binomial')
p3a_c.bin=glm(pres.approve.w2.bin~topic.w2+president.w2+condition.w2+pres.approve.w1+sex+age+income+ethnicity+copartisanship*outcome,data=dt,family='binomial')

outcome_comparisons_party=function(model,data,simulations=1000,seed=6658) {
  set.seed(seed)
  betas=slot(sim(model,simulations),'coef')
  covariates=names(model$model)
  results=matrix(data=NA, ncol=7,nrow=3)
  
  # add columns for factor coefs
  D=data[complete.cases(data[,covariates]),]
  D=cbind(D,model.matrix(~topic.w2-1,data=D))
  D=cbind(D,model.matrix(~president.w2-1,data=D))
  if ('copartisanship' %in% covariates) {D=cbind(D,model.matrix(~copartisanship-1,data=D))}
  if ('topic.approve.w1' %in% covariates) {D=cbind(D,model.matrix(~topic.approve.w1-1,data=D))}
  if ('pres.approve.w1' %in% covariates) {D=cbind(D,model.matrix(~pres.approve.w1-1,data=D))}
  if ('manip.done.w1' %in% covariates) {D=cbind(D,model.matrix(~manip.done.w1-1,data=D))}
  if ('copartisanship' %in% covariates) {D=cbind(D,model.matrix(~partisanship-1,data=D))}
  if ('sex' %in% covariates) {D=cbind(D,model.matrix(~sex-1,data=D))}
  if ('income' %in% covariates) {D=cbind(D,model.matrix(~income-1,data=D))}
  if ('ethnicity' %in% covariates) {D=cbind(D,model.matrix(~ethnicity-1,data=D))}
  if ('education' %in% covariates) {D=cbind(D,model.matrix(~education-1,data=D))}
  D=cbind(D,model.matrix(~condition.w2-1,data=D))
  
  # intercept + covariates 
  controls=cbind(1,D$`topic.w2contractor pay`,D$`topic.w2farm payments`,D$`topic.w2foreign aid and abortion`,D$`topic.w2foreign workers`,D$`topic.w2gun violence research`,D$topic.w2LGBT,D$`topic.w2policing and military surplus`,D$`topic.w2public lands`,D$`topic.w2Russian sanctions`,D$`topic.w2student loans`,D$topic.w2trade,D$`topic.w2water rules`,D$`topic.w2weapon sales to Saudi Arabia`,D$topic.w2wildlife,D$president.w2Trump,D$condition.w2congress,D$condition.w2order)
  if ('topic.approve.w1' %in% covariates) {controls=cbind(controls,D$topic.approve.w1Disapprove,D$`topic.approve.w1Somewhat disapprove`,D$`topic.approve.w1Neither approve or disapprove`,D$`topic.approve.w1Somewhat approve`,D$`topic.approve.w1Approve`,D$`topic.approve.w1Strongly approve`)}
  if ('pres.approve.w1' %in% covariates) {controls=cbind(controls,D$pres.approve.w1Disapprove,D$`pres.approve.w1Somewhat disapprove`,D$`pres.approve.w1Neither approve or disapprove`,D$`pres.approve.w1Somewhat approve`,D$`pres.approve.w1Approve`,D$`pres.approve.w1Strongly approve`)}
  if ('manip.done.w1' %in% covariates) {controls=cbind(controls,D$manip.done.w1Disagree,D$`manip.done.w1Somewhat disagree`,D$`manip.done.w1Neither agree nor disagree`,D$`manip.done.w1Somewhat agree`,D$`manip.done.w1Agree`,D$`manip.done.w1Strongly agree`)}
  if ('sameparty.w1' %in% covariates) {controls=cbind(controls,D$sameparty.w1)}
  if ('sex' %in% covariates) {controls=cbind(controls,D$sexFemale)}
  if ('age' %in% covariates) {controls=cbind(controls,D$age)}
  if ('income' %in% covariates) {controls=cbind(controls,D$`income2: $25,000 to $50,000`,D$`income3: $50,001 to $75,000`,D$`income4: $75,001 to $100,000`,D$`income5: More than $100,001`)}
  if ('ethnicity' %in% covariates) {controls=cbind(controls,D$ethnicityBlack,D$`ethnicityNative American`,D$`ethnicityOther/Decline to State`,D$ethnicityWhite)}
  if ('education' %in% covariates) {controls=cbind(controls,D$`educationSome College or Vocational`,D$`educationB.A. or B.S.`,D$`educationPost-graduate or Higher`)}
  
  # new data to simulate over
  d.success.opp=cbind(controls,0,1,0,0,0)
  d.fail.opp=cbind(controls,0,1,1,0,1)
  
  d.success.cop=cbind(controls,1,0,0,0,0)
  d.fail.cop=cbind(controls,1,0,1,1,0)
  
  d.success.ind=cbind(controls,0,0,0,0,0)
  d.fail.ind=cbind(controls,0,0,1,0,0)
  
  # simulations
  p.win.opp=apply((1/(1+exp(-as.matrix(d.success.opp)%*% t(betas)))),1,as.vector) 
  p_win.opp=apply(p.win.opp,1,mean)
  p.fail.opp=apply((1/(1+exp(-as.matrix(d.fail.opp)%*% t(betas)))),1,as.vector) 
  p_fail.opp=apply(p.fail.opp,1,mean)
  
  p.win.cop=apply((1/(1+exp(-as.matrix(d.success.cop)%*% t(betas)))),1,as.vector) 
  p_win.cop=apply(p.win.cop,1,mean)
  p.fail.cop=apply((1/(1+exp(-as.matrix(d.fail.cop)%*% t(betas)))),1,as.vector) 
  p_fail.cop=apply(p.fail.cop,1,mean)
  
  p.win.ind=apply((1/(1+exp(-as.matrix(d.success.ind)%*% t(betas)))),1,as.vector) 
  p_win.ind=apply(p.win.ind,1,mean)
  p.fail.ind=apply((1/(1+exp(-as.matrix(d.fail.ind)%*% t(betas)))),1,as.vector) 
  p_fail.ind=apply(p.fail.ind,1,mean)
  
  # store
  results[1,1]<- quantile(p_fail.opp-p_win.opp,.0083)
  results[1,2]<- quantile(p_fail.opp-p_win.opp,.025)
  results[1,3]<- mean(p_fail.opp-p_win.opp)
  results[1,4]<- quantile(p_fail.opp-p_win.opp,.975)
  results[1,5]<- quantile(p_fail.opp-p_win.opp,.9916)
  results[1,6]<- 'Opposition'
  results[1,7]<- covariates[1]
  
  results[2,1]<- quantile(p_fail.cop-p_win.cop,.0083)
  results[2,2]<- quantile(p_fail.cop-p_win.cop,.025)
  results[2,3]<- mean(p_fail.cop-p_win.cop)
  results[2,4]<- quantile(p_fail.cop-p_win.cop,.975)
  results[2,5]<- quantile(p_fail.cop-p_win.cop,.9916)
  results[2,6]<- 'Copartisans'
  results[2,7]<- covariates[1]
  
  results[3,1]<- quantile(p_fail.ind-p_win.ind,.0083)
  results[3,2]<- quantile(p_fail.ind-p_win.ind,.025)
  results[3,3]<- mean(p_fail.ind-p_win.ind)
  results[3,4]<- quantile(p_fail.ind-p_win.ind,.975)
  results[3,5]<- quantile(p_fail.ind-p_win.ind,.9916)
  results[3,6]<- 'Independents'
  results[3,7]<- covariates[1]
  
  return(results)
}

other_outcomes=data.frame(lb.MC=numeric(9),lb.95=numeric(9),me=numeric(9),ub.95=numeric(9),ub.MC=numeric(9),grp=character(9),dv=character(9))
other_outcomes[1:3,1:7]=outcome_comparisons_party(t3a_c.bin,dt)
other_outcomes[4:6,1:7]=outcome_comparisons_party(p3a_c.bin,dt)
other_outcomes[7:9,1:7]=outcome_comparisons_party(v3a_c,dt)

# manip for figure
final_outcomes=rbind(other_outcomes,outcomes)

final_outcomes$lb.MC=as.numeric(final_outcomes$lb.MC)
final_outcomes$lb.95=as.numeric(final_outcomes$lb.95)
final_outcomes$me=as.numeric(final_outcomes$me)
final_outcomes$ub.95=as.numeric(final_outcomes$ub.95)
final_outcomes$ub.MC=as.numeric(final_outcomes$ub.MC)

final_outcomes$dv=as.factor(final_outcomes$dv)
levels(final_outcomes$dv)[levels(final_outcomes$dv)=='sameparty.w2'] <- 'Same-Party\nVote'
levels(final_outcomes$dv)[levels(final_outcomes$dv)=='topic.approve.w2.bin'] <- 'Topic\nApproval'
levels(final_outcomes$dv)[levels(final_outcomes$dv)=='pres.approve.w2.bin'] <- 'Job\nApproval'

f2_b = ggplot(final_outcomes) +
  geom_hline(yintercept = 0,
             colour = gray(1 / 2),
             lty = 2) +
  geom_linerange(aes(x = dv, ymin = lb.MC, ymax = ub.MC),
                 lwd = 1.1,
                 position = position_dodge(width = 1 / 2)) +
  geom_pointrange(
    aes(
      x = factor(
        dv,
        levels = c('Same-Party\nVote', "Job\nApproval", "Topic\nApproval")
      ),
      y = me,
      ymin = lb.95,
      ymax = ub.95
    ),
    lwd = 2,
    size = 2.5,
    position = position_dodge(width = 1 / 2),
    shape = 21,
    fill = "WHITE"
  ) +
  geom_text(
    aes(
      x = factor(
        dv,
        levels = c('Same-Party\nVote', "Job\nApproval", "Topic\nApproval")
      ),
      y = me,
      label = round(me, 2)
    ),
    size = 3,
    position = position_dodge(width = 1 / 2)
  ) +
  facet_wrap( ~ grp, ncol = 4, scales = 'free_x') +
  coord_flip() + theme_bw() + xlab('') + ylab('') +
  theme(
    text = element_text(family = 'Palatino'),
    legend.position = 'bottom',
    legend.title = element_text(size = 9),
    strip.background = element_rect(
      color = "black",
      fill = "#ffffff",
      linewidth = 1.5,
      linetype = "blank"
    )
  ) +
  guides(color = guide_legend(
    override.aes = list(shape = 19, size = 1),
    reverse = T,
    title.position = "right"
  ))

ggsave("figures-tables/sameparty_f3.pdf",
       f2_b,
       width = 7.2,
       height = 4,
       units = 'in')


